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We study the dynamics of weakly coupled non-abelian plasmas within the frameworks of classical- 
statistical lattice gauge-theory and kinetic theory. We focus on a class of systems which are highly 
occupied, isotropic at all times and initially characterized by a single momentum scale. These 
represent an idealized version of the situation in relativistic heavy ion-collisions in the color-glass 
condensate picture, where on a time scale 1/Qs after the collision of heavy nuclei a longitudinally 
expanding plasma characterized by the saturation scale Q s is formed. Our results indicate that the 
system evolves according to a turbulent Kolmogorov cascade in the classical regime. Taking this 
into account, the kinetic description is able to reproduce characteristic features of the evolution 
correctly. 
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I. INTRODUCTION 

Understanding the mechanisms underlying the fast 
thermalization observed in relativistic heavy-ion colli- 
sions [TJ provides one of the biggest theoretical challenges 
in our current understanding of the experiments per- 
formed at RHIC and the LHC g]. Though there has 
been significant progress in the development of ab-initio 
calculations [3J, these have not yet been able to address 
the problems of thermalization and isotropization [3] . On 
the other hand several authors have studied thermaliza- 
tion of non-abelian plasmas in a more general context at 
strong [Jj and weak coupling [5HS]. We follow the same 
approach and study the problem of thermalization from 
a more general point of view, while focussing on a class 
of systems which share important features with the non- 
equilibrium stage of relativistic heavy ion collisions. We 
will work at weak coupling a <C 1 and consider a class 
of homogenous, isotropic sysyems which initially have a 
parametrically large occupation 



/(P) 



for |p| < Q , /(|p| > Q) « 1, (1) 



where < c < 1, and are characterized by a single di- 
mensionful scale Q. As discussed in Ref. [5J this mimics 
the behavior in heavy-ion collisions at times t ~ Q^ 1 
after the collision of heavy nuclei. The main difference 
is the longitudinal expansion of the system, which we 
will neglect in the following. This greatly simplifies the 
problem, as the system can be considered isotropic at all 
times, whereas in the longitudinally expanding case the 
system remains anisotropic on large time scales [5J [S]. 
Thermalization of such a system has been studied in a 
kinetic theory framework in Refs. E] and first results 
from classical-statistical lattice simulations have been 
presented in Ref. [JO]. While the authors of Refs. [7] 
and [5] agree on the way thermalization proceeds in these 
type of overoccupied systems, numerical simulations have 
raised the question whether alternatively thermalization 
may proceed as a turbulent cascade [JO] [JJJ . The latter 
are characterized by non-thermal fixed points of classi- 
cal field theories, which can be associated to stationary 
transport of conserved quantities [JJ] [T3J. The occur- 
rence of these solutions is well known for a variety of 



systems including early universe cosmology |14j and cold 
atomic gas systems [T5T - H7] . and turbulent behavior has 
also been discussed in the context of non-abelian gauge 
theories [TOl [TT1 ITM2"0] . In this paper we will focus on the 
dynamical evolution of the system and refer to the liter- 
ature for analytic studies of turbulence [TOHTSl [T8ll20| . 
Our strategy is to combine kinetic theory and classical 
statistical lattice simulations in order to provide a simple 
understanding based on parametric estimates while hav- 
ing control over non-perturbative and non-equilibrium ef- 
fects. We start with a kinetic theory analysis in Sec. |TIJ 
where we briefly review the discussion in Refs. [7] [5] and 
extend it to account for the presence of a non-thermal 
fixed point. In Sec. Ill we present results from classical- 



statistical lattice simulations, which indicate the occur- 
rence of a turblent cascade with exponent k ~ 4/3 at late 
times, as observed in Ref. We then summarize our 

results and conclude with Sec. HVl 



II. EVOLUTION IN KINETIC THEORY 

The kinetic evolution of systems characterized with an 
initial distribution as in Eq. has been studied in Refs. 
[7J [8] and we will adopt the notation of Ref. [5] to an- 
alyze the evolution of the system. To describe the evo- 
lution, one has to take into account the effects of elastic 
and inelastic scattering, as well as the interaction of hard 
and soft excitations. We start with a discussion of elas- 
tic scattering, which in this case is dominated by scat- 
tering of hard particles with small momentum transfer 
|5J . In the small angle approximation, where the aver- 
age momentum of hard excitations is much bigger than 
the average exchanged momentum, the process appears 
as a 'random walk' in momentum space, controlled by 
the momentum diffusion parameter (c.f. Ref. [7J [5]) 



d 3 p 
(2^)3 



/(p)[l+/(p)] 



(2) 



Similarly one can estimate the effects of inelastic pro- 
cesses and interactions of hard particles with soft gauge 
field modes below the Debye scale. The important ob- 
servation is that for overoccupied, isotropic systems all 
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processes are parametrically as efficient as elastic scat- 
tering [Jj 18] such that 



q so ft 



(3) 



This greatly simplifies the discussion of the kinetic evo- 
lution of the system, since we do not have to distinguish 
different regimes where either of the processes dominate. 
FromEq. ^ we can also obtain the rate at which a par- 
ticle with momentum p exhibits a momentum transfer of 
the same order. This rate is parametrically given by 



r ei (p) 



<lel_ 
r>2 



(4) 



such that soft particles experience large angle scatterings 
at a higher rate as compared to hard particles. Given the 
enhanced scattering rates of soft excitations, we expect 
that the soft tail of the distribution can always 'equili- 
brate' to a distribution which is dictated by the dynamics 
of the hard modes. We will assume in the following that 
this distribution is described by a power law, where 



/(P) 



for p < A 



(5) 



up to a momentum cutoff A above which occupancies are 
negligible. Here the exponent c characterizes paramet- 
rically the occupation of the system, and we will only 
consider the case < c < 1 of initially over-occupied 
systems. The exponent k characterizes the shape of the 
distribution, while k — 1 corresponds to a thermal shape, 
the values k = 3/2 (c.f. Ref. [TO]) and k = 4/3,5/3 (c.f. 
Ref. [llj ) are known to appear for systems which exhibit 
Kolmogorov wave turbulence [j] The scales A and A s in 
Eq. ([5]) both depend on time and determine the evolution 
of the system. At initial time to ~ Q _1 the two scales 
coincide and one finds A ~ A s ~ Q. In order to extract 
the time evolution of the scales A and A s we first note 
that for k < 3/2 the momentum diffusion parameter is 
parametrically given by 



,2-2c 



A 



(6) 



In contrast if one considers k > 3/2 the integral in Eq. 
([2]) is dominated by infrared contributions indicating a 
break down of the kinetic description. We will therefore 
limit our discussion to the case n < 3/2, where we ex- 
pect a kinetic description to apply. As we will see shortly 
in Sec. |III[ this situation is also realized in our lattice 
simulations. From Eq. Q we can infer the time scale 
^start on which the momentum of initially hard excita- 
tions changes appreciably as 



Q 



(7) 



1 The appearance of the exponent k = 3/2 can be explained with 
the existance of a classical background field I10II13| . This aspect 
is discussed in more detail in Ref. |10| in the context of the recent 
debate on condensate formation out-of equilibrium 8 21. 



Before the time i s tart the system will develop its soft tail, 
which then moves outwards to higher momenta until at 
times t ~ t s tart changes in the distribution of hard exci- 
tations start to take place. In the regime t ^ ^start one 
can find a self-consistent scaling solution by requiring the 
evolution of the hard scale to be governed by momentum 
diffusion such that (c.f. Ref. [3 H]) 



dt 



-A 7 



a 



A' 5 



(8) 



We note that for non-expanding systems the ratio of soft 
and hard scale (A S /A) K is fixed by energy conservation 



e ~ a 



such that the relation 



A 4 ~cV 



c n 4 



(9) 



(10) 



holds at all times of the evolution. Combining the above 
conditions then yields the evolution of the scales A and 
A s for times t > i s tart to be 



A ~ Q 



t 



1/7 



, A S ~Q 



t 



(i-4/«)/7 



(11) 



which agrees with the results obtained in Refs. [7J |S] 
for the case K = 1 of a (quasi-)thermal distribution. We 
note that the evolution of the hard scale A as well as 
the evolution of the occupancies of hard modes (A S /A) K 
turn out to be independent of the exponent k. This 
behavior is a direct consequence of energy conservation. 

The evolution of the scales A and A s proceeds in 
this way until at some points the occupancies of hard 
modes become 0(1), where quantum effects become 
important. These will have the effect of driving the 
system towards its unique thermal fixed point, and 
hence complete the thermalization process. We therefore 
expect thermalization of the system to occur as two step 
process, where the first regime is characterized by the 
above scaling solutions. The timescale t c i iangc for enter- 
ing the second regime, where the 'rest-thermalization' 
takes place is parametrically given by 



change 



Q 



"l a -2+c/4 



(12) 



which is independent of the value of the exponent k. 
However for k / 1 even at times t ~ Change the system 
still deviates significantly from its equilbirium distribu- 
tion and one has to take into account also the time for the 
'rest-thermalization' to occur. If the latter is assumed to 
be parametrically faster than our previous estimate one 
obtains t cq ~ t c hangc as the final result. In this situation 
one recovers the estimate of Refs. [7J U| for the overall 
thermalization time. 
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III. EVOLUTION IN CLASSICAL-STATISTICAL 
LATTICE GAUGE-THEORY 

In this section we present the results from classical 
statistical lattice simulations of SU(2) Yang-Mills 
theory. The primary goals are to analyze whether 
the above scaling solutions are realized in a genuine 
non-equilibrium approach and to determine the scaling 
exponent k. 

When comparing the lattice results to the kinetic 
theory evolution, we have to overcome the challenge of 
identifying the quasi-particle excitations discussed in the 
kinetic theory framework. This is not unambiguous since 
the notion of quasi-particles is a perturbative concept 
and there is no general correspondence to two-point 
correlation functions. However one can typically exploit 
the perturbative behavior of the two-point correla- 
tion functions to establish a quasi particle picture, 
in situations where this is applicable. With regard 
to non-abelian gauge-theories the problem becomes 
even more severe, due to the problems associated with 
gauge-invariance. Typically the particle number defined 
from two-point correlation functions will not be a 
gauge invariant quantity and it is therefore important 
to consider a gauge which has a clear interpretation 
in terms of physical degrees of freedom. In practice 
the choice of gauge is limited by the necessity to use 
temporal axial gauge (A t = 0) in the classical-statistical 
lattice simulations. However the gauge fixing A t — 
is incomplete such that one has the residual freedom 
to perform time-independent gauge transformations. 
Fixing the residual gauge by an additional constraint 
eliminates gauge-fluctuations and can be used to define 
quasi-particle observables. We follow Ref. [TU] and 
address this problem by implementing the Coulomb type 
constraint V A — 0, whenever we are interested in gauge- 
dependent observables. We then study the evolution of 
the electric and magnetic two-point correlation functions 



(BB)(t,p)= I 
(EE)(t,p)= f 

■' x-y 



(Bf(f,x)Bf(t,y)) 
(JV?-l)(d-l) 
x)iff(t,y)) 
(JV?-l)(d-l) ' 



p(x-y) 



-ip(x-y) 



(13) 



(14) 



here written in terms of continuum variables, where N c 
is the number of colors, d = 3 is the number of spa- 
tial dimensions and (.) denotes classical-statistical aver- 
aging^] The classical-statistical lattice implementation is 
describe in detail in Ref. [22] and we use the so-called 
'Los Alamos' gauge-fixing procedure [23]- The initial 
conditions in our simulations are chosen to mimic the 



2 The normalization factor (N^ — l)(d — 1) corresponds to the 
number of massless physical gluons. This normalization implic- 
itly assumes that the unphysical degrees of freedom are removed 
by the gauge-fixing procedure. 
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FIG. 1. (color online) Spectrum of excitations obtained from 
the correlation functions, (EE)/\p\ and {BB)/\p\ at different 
times of the evolution. One observes a clear power-law de- 
pendence for the correlation functions (EE) and (BB) , which 
extends over approximately one decade. As time proceeds the 
tail of the power law propagates towards the ultraviolet, while 
the amplitude of the distribution decreases. 



quasi-particle picture in Eq. ([I]) and described in more 
detail in the appendix. The simulations are performed on 
N = 64,96,128 hypercubic lattices with three different 
values of the lattice spacing Qa — 0.33,0.66, 1.0 to gain 
systematic control over lattice spacing and finite volume 
effects. If not stated otherwise results are presented for 
N = 96 and Qa = 0.66 lattices. 



1. Highly overoccupied systems 

r/o( P )~a- 1 ; 

We first study the time evolution of the correlation 
functions defined in Eq. (13) and (14). The results are 



presented in Fig. [T] where we show the spectrum of the 
correlation functions (EE)(t, p)/|p| and (BB)(t, p)/|p| 
at different times of the evolution. The factor l/|p| is 
chosen to produce dimensionless quantities, which can 
be interpreted as an occupation numbers. From Fig. [T] 
one observes the build up of a power law distribution, 
which subsequently decreases in amplitude while slowly 
moving out towards higher momenta. We observe that 
towards later times the electric and magnetic correlation 
functions (EE) /\p\ and (BB)/\p\ show a consistent scal- 
ing behavior over a wide range of momenta. In order to 
establish a more direct comparison with the kinetic the- 
ory discussion we define the effective occupation number 
from the electric field correlation function as 



f(t,p) 



(EE)(t,p) 



(15) 
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FIG. 2. (color online) Scaling exponent k extracted from the 
spectrum of the correlation function f(t,p) ~ {EE)/\p\ at 
different times of the evolution. The different symbols corre- 
spond to the two different extraction methods. The red filled 
squares correspond to a global fit in the momentum range 
p — 0.3 - 1.0 Q, whereas the black empty squares correspond 
to average of local fits as described in the text. Over a large 
time-scale the exponent is consistent with the Kolmogorov 
turbulence exponent k = 4/3 [11] , 



which we will use in the following to further analyze the 
evolution. 

While the power law dependence of the distribu- 
tion function f(t,p) can already be observed from Fig. 
[I] we are interested also in the exponent k of the power 
law. This is important to distinguish the cases where 
the observed spectrum corresponds to a quasi-thermal 
evolution (k = 1) or a turbulent Kolmogorov cascade 
(k = 4/3,5/3; 3/2). To extract the exponent from the 
data we perform a series of least-square fits at each time 
slice, with a distribution function of the form 



f(t,p) 



a 



A s (t) 



K(t) 



(16) 



This procedure yields the two parameters n(t) and 
(A s / Q) K (t) which are shown in Fig. [2] and [3] respectively. 
To estimate the error of this procedure we perform the 
analysis in two different ways: In the first case we simply 
consider the result of a global fit along with its error in 
the momentum range p/Q — 0.3 — 1.0, where scaling is 
observed. The results of this procedure are shown as 
red squares in Fig. [2] In the second case we divide the 
data in momentum bins and perform a separate fit for 
each bin. We then extract the average exponent and its 
error in the scaling region. The results correspond to the 
black points in Fig. [2] One observes that the results for 
k of the two procedures agree, while the second method 
provides a more reliable estimate of the error. We find 
that at early times the spectrum features a power law 
exponent of k ~ 3/2, which is consistent with previous 
observations |10) . At later times we observe a hardening 
of the spectrum where the data is in favor of k ~ 4/3 



FIG. 3. (color online) Time evolution of the soft scale 
(A S /Q) K extracted from fits of the spectra at different times 
(red squares) and from the time evolution of the occupation 
number for modes with p ~ Q (black squares). The blue 
dashed line corresponds to a power-law behavior (Qi)' K_4 " )/ ' 7 
as expected from the kinetic theory analysis, where we as- 
sumed k = 4/3 as a constant in time. 



over a large time-scale. The appearance of the exponent 
k = 3/2, may be attributed to a transient condensation 
phenomenon, which does not persist on parametrically 
large time scales [TUl 121] • The thermal value n = 1 is 
clearly ruled out by the data over the entire simulation 
time. 

As discussed in Sec. |n] in the kinetic theory 
framework, the appearance of a non-thermal exponent 
(k 7^ 1) has an immediate impact on the evolution of 
the soft scale A s , which is expected to display a slower 
evolution for larger values of k. To investigate whether 
this is supported by the lattice data, we extract the time 
evolution of the soft scale (A s /Q) K (t) by two different 
procedures. The results are on display in Fig. [3] In the 
first case we use the results of the previous fit procedure 
to directly obtain the quantity (A S /Q) K as a function of 
time. This corresponds to red squares in Fig. [3] In the 
second case we investigate the evolution of modes with 
|p| ~ Q, where f(t, |p| ~ Q) ~ a~ c (A s /Q) K according 
to Eq. (16). This is shown by black squares in Fig. [3] 



The results from both methods agree and one observes 
a clear power law dependence. The blue dashed line 
corresponds to the kinetic theory estimate (c.f. Sec. |n| 



A s (t) 
Q 



(Qt) 



(N-4)/7 



(17) 



where we assumed k = 4/3 to be independent of time. 
One observes good overall agreement with the data. We 
also performed fits to the extract the scaling exponent 
from the time evolution of (A S /Q) K . Using the same 
methods as introduced previously and assuming the 
relation (17) this yields n ~ 1.3 ±0.1, which is consistent 
ues shown in Fig. [2j 



with the va 
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FIG. 4. (color online) Rescaled moments of the distribution 
function as a function of the rescaled momentum variable. 
The different colors and symbols correspond to different evo- 
lution times. The scale factor is s = t/to, where we chose 
Qto = 1000 as normalization. The scaling exponents are sen- 
sitive to the evolution of the hard scale A(t) and we used the 
values from the kinetic theory analysis in Sec. [II] The results 
are obtained on N = 96 lattices with Qa — 0.33. 



So far we have only analyzed the behavior of the 
spectrum in a momentum region where it shows a clear 
power law dependence. However one of the key features 
of the evolution is the propagation towards higher 
momenta, which ultimately leads to thermalization of 
the system. As discussed in Sec. [TT] the kinetic theory 
framework provides a strong prediction of how this 
evolution proceeds, which is insensitive also to the power 
law exponent k. To extract the evolution of the cut-off 
scale A(t) from the lattice data, we exploit the fact 
that the spectrum shown in Fig. [T] shows a self-similar 
behavior. To illustrate this we first note that we can 



parametrize the entire spectrum as 
'A s (i)V 



f(t,p) 



C(t, |p|/A(f)), (18) 



where the first part corresponds to the infrared power-law 
and the second part regulates the ultra-violet behavior 
such that C(t, |p|/A(t)) = 1 for |p| < A(t), whereas it 
drops off quickly for |p| > A(t) as observed in Fig. [I] 
Here k ~ 4/3 is assumed to be constant in time. If we 
assume further that the shape of the cut-off function is 
independent of time, i.e. 



C(t,|p|/A(t))^C(|p|/A(t)) 



(19) 



and the scales A(t) and A s (t) evolve according to a power 
law in time, i.e. 



A(t) 
A,(t) 



Q (tAstart)" 
Q (iAstart)' 3 



(20) 
(21) 



we find that the distribution function is self-similar in the 
sense that for s > one finds 



f(st,s a p) = s 



f(t,p) ■ 



(22) 



Here it is crucial to rescale momenta according to s a , 
whereas the evolution time scales with s in order to 
reproduce the correct cut-off behavior. Therefore Eq. 
) is particularly sensitive to the scaling of the hard 
e A(t) an can be used to extract the scaling exponent 



(22 



sea 
a. 



We note that the assumption ( 19 1 only concerns the 



behavior of p > A(t), while for p < A(t) this property 
is satisfied by construction. Hence the condition (19) 



reflects the assumption, that the fall-off above the cut-off 
scale is universal and we will turn back to this point 
when we contrast the above analysis with lattice data. 
Also the scaling behavior of A s (t) readily emerges from 
Fig. |3j where we confirmed (3k = (k — 4)/7 in agreement 
with the analysis in Sec. |Hj where we obtained the 
exponents a = 1/7 and /? = (1 — 4/k)/7. 

In order to investigate to which degree this self- 
similarity is featured by the lattice data, we invert Eq. 
( 22 ) such that we expect the spectra at different times 
to collapse on a universal curve. As we are particularly 
interested in the evolution of the hard scale A(i), we find 
it convenient to consider moments of the distribution 
function 



/W(t,p) = |pp/(i,p) 



(23) 



which are more sensitive to hard modes than the distribu- 
tion function itself. Setting s = t/to the scaling relation 
for the moments f( n \t,p) then takes the form 

/ (n) (p) = (iAo) (4 - n)/ 7 (n) (i, (t/t ) 1/7 p) (24) 

where we used the values of a, (3 from Sec. [H] and ab- 
breviated /o n) (p) = / (n) ( i o,p) for arbitrary t > t stalt - 
The rescaled moments of the lattice data are shown in 
Fig. [4] for n = 2,3 and 4 and times Qt = 1000 - 8000, 
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Momentum p/Q 



FIG. 5. (color online) Spectrum of excitations for no = 0.05 
at different times of the evolution. In addition to A 7 = 96 
and Qa — 0.66 data (open symbols) we also show results with 
reduced statistics for N = 128 and Qa — 1 (filled symbols). 
One observes that the evolution of modes with p ~ Q is much 
slower as compared to the no = 1 case. As a consequence 
only the soft sector is subject to changes due to interactions 
at early times, while at later times hard modes are also af- 
fected. In the transition region the spectrum shows a thermal 
shape k ~ 1 at some point of the evolution. Nevertheless the 
evolution at late times proceeds via a turbulent cascade with 
k ~ 4/3. This behavior is indicated by the black dashed lines. 



FIG. 6. (color online) Evolution of the soft scale (A S /Q) K 
normalized by the initial occupation no as a function of 
the normalized time t/tstart for different initial occupations 
no = 0.05,0.1,0.2,0.5 and 1 (bottom to top). The quan- 
tity (A S /Q) K is obtained from the occupation of modes with 
p ~ 0.9Q. The black dashed line portraits the result of the ki- 
netic theory analysis (c.f. Sec. [TTp . At late times t ^> fgtart one 
observes a power law dependence independent of the choice 
of no in the considered range of parameters. 



for all values of no considered here. 



where we used Qto = 1000 as the normalization. As 
the moments are more sensitive to the hard tail of the 
distribution, we used N — 96 and Qa = 0.33 in our 
simulations to achieve a larger momentum cut-off . One 
observes that the rescaled data collapses approximately 
on a single curve. This provides strong evidence that the 
hard scale A(t) indeed follows a power law behavior with 
A(t) ~ Q (Qt) 1 / 7 as discussed in Sec. |nj In particular 
the position and amplitude of the peak, which are very 
sensitive to a — 1/7, coincide for all curves. Above the 
peak we observe minor deviations from the universal be- 
havior. This can be attributed to a non-universal fall-off 
of the cutoff function for p > A(t), however we can also 
not completely rule out the presence of lattice artifacts 
in the high momentum region. 



2. Overoccupied systems 

C/ (p)~Q- C , 0<C<1) 

We now extent our previous analysis to systems with 
initial occupancies which are still parametrically large 
but now smaller than a -1 . We will in the following 



denote the initial occupancy by no 



a 



and we study 



systems with n — 0.05,0.1,0.2,0.5 for a fixed value of 
the strong coupling constant a = 10~ 6 such that the 
classicality condition no a is always satisfied. We will 
first discuss the case no = 0.05 as an example and then 
turn to a global analysis of all data. We note that the 
qualitative behavior on large time scales is very similar 



We proceed as previously and first study the evo- 
lution of the spectrum of the correlation functions. In 
Fig. [5] we present snapshots of the spectrum of f(t, p) as 
defined in Eq. (15) for initial occupation no = 0.05. At 



early times we observe an increase of occupancies in the 
soft sector, whereas the occupation for modes with p ~ Q 
remains almost unaffected. This is in accordance with 
our expectations from the kinetic theory analysis (c.f. 
Sec. [n]). The evolution at late times proceeds, similarly 
to the case no = 1, as a turbulent cascade (k ~ 4/3) 
towards the ultraviolet. The transition region between 
early and late times is characterized by a softening of the 
spectrum, where the initially flat distribution changes 
its shape to a power-law distribution with exponent 
k ~ 4/3. Even though the classical thermal value 
k ~ 1 is featured at some time during this transition, 
the system subsequently continues to evolve towards 
k ~ 4/3 at later times. This provides yet another strong 
piece of evidence that a non-thermal exponent k ~ 1 is 
indeed favored during the far from-equilibrium evolution. 

When comparing the results in Fig. [5] with the 
evolution in the case no = 1, we observe that the 
evolution proceeds significantly slower in the case of 
reduced occupancy. This is expected given the smaller 
interacting rates for lower occupancies. The kinetic 
theory analysis presented in Sec. [TT] here makes a 
strong statement about the occupancy dependence, 
based on the parametric dependence of the momentum 
diffusion parameter q ~ n^Q 3 , which controls the rate of 
interactions at the hard scale. Accordingly the natural 
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time scale for the system to evolve is therefore no 
longer set by Q _1 but rather by the 'scattering time' 
igtart ~ Q ln o 2 which initially controls the rate of 
interactions at the hard scale (see also Ref. [?.). In 
order to analyze whether this behavior is displayed 
by the lattice data, we chose to investigate the time 
evolution of the soft scale (A S /Q) K for different initial 
occupancies. This quantity can be obtained from the 
occupation of modes with p ~ 0.9 Q and is shown in 
Fig. [6] as a function of time in units of the 'scattering 
time'. One observes that around t ~ istart the transition 
from an approximately constant behavior to a power-law 
decay takes place. The time dependence at late times 
t S> tstart is well described by a power law, where 
(A S /Q) K - 0Astart) (K - 4)/7 with k ~ 4/3 as indicated 
by the black dashed line in Fig. [6| Concerning the 
dependence on the initial occupancy no, we find that the 
expected scaling holds approximately as the data shown 
in Fig. [6] nearly collapses on a single curve. However 
one does observe a residual dependence of the overall 
amplitude of (A S /Q) K on the initial occupancy, which we 
found to be rather sensitive to the details of the initial 
conditions. 



IV. SUMMARY AND DISCUSSION 

In this paper we studied the evolution of weakly 
coupled non-abelian plasmas, which are characterized 
by parametrically large occupancies. In Sec. [TT] we 
extended the kinetic theory discussion of Refs. [7J |S] 
to account for the presence of a non-thermal fixed 
point, characterized by a power law distribution with 
exponent n. We found that for k < 3/2 the dynamics 
is dominated by hard modes and the analysis in Refs. 
|H] can be extended to non-thermal systems in a 
straightforward way. In particular the evolution of the 
hard scale A and the occupancies at the hard scale 
(A S /A) K , which are relevant to estimate the thermaliza- 
tion time are unaffected by this generalization. In Sec. 
III| we studied the time evolution of the system within 
classical-statistical real-time lattice simulations. We 
presented strong evidence that the evolution proceeds 
as a turbulent cascade with exponent k ~ 4/3, where 
the system exhibits a self-similar behavior. This is 
further supported by the results reported in Ref. [TT] . 
where k = 4/3 has been obtained from both numerical 
and analytical considerations. When taking this into 
account in the kinetic description, we found good overall 
agreement with the lattice data. In summary our results 
suggest a 'turbulent thermalization' mechanism, which is 
very similar to observations in early-universe cosmology 
[14] and systems of cold-atomic gases [ISHIZ]- 



In view of the early time dynamics of relativistic 
heavy-ion collisions it is important to take into account 
the longitudinal expansion of the system. The latter 
renders the system anisotropic on large time scales 



and additional effects such as plasma instabilities are 
expected to play an important role [9]. It is therefore 
not possible to directly transfer our results to expanding 
systems and separate numerical studies will have to 
be performed. While significant progress has recently 
been made in simulations of expanding systems in 
scalar quantum field theories |24j . present simulations 
in non-abelian gauge theories 01 have not yet been able 
to resolve this question and we expect future studies to 
clearify the situation. 

Acknowledgement: The author likes to thank J. 
Berges, L. McLerran, D. Sexty and in particular R. 
Vcnugopalan for insightful discussion and collaboration 
on related projects. The author also thanks Brookhaven 
National Lab for hospitality during his stay, where this 
work was performed. This work was supported in part 
by BMBF grant 06DA9018 and by HGS-HIRe for FAIR. 



APPENDIX A: INITIAL CONDITIONS 

The initial conditions for our classical-statistical lattice 
simulation are chosen to mimic a quasi particle picture, 
where we consider a superposition of modes CJA,a(iojP) 
labeled by the color index a = 1,...,N^ — 1, the po- 
larization A = 1,2 and spatial momentum p. The 
mode functions a\, a (to,p) are the ones of the free the- 
ory, which satisfy the gauge-condition p i a\ a (to ) p) = 
as well as the abelian part of the gauss-law constraint 
Pi@t o,\(t, p)| t=to = individually for each mode. The 
initial occupation of modes is then determined by the 
relations 



(a x , a (t,p)a* x , b {t',p))\ t=t , =tg 
dt(a\, a (t,p)a* x , jb (t'p))\ t=t , =tg 
d t dt> (o A , (t, p)a* x , jb {t'p)) \ t=t , =t 



SabS\\'N p /2u} p 


6 ab 5\yuj p N p /2 

(25) 



where N p = fo (p) + 1/2 is the initial distribution of ex- 
citations. We consider systems with parametrically large 
occupation according to 



/o(p) = 



(26) 



where 9{x) is the Heavyside step function, a <C 1 is the 
strong coupling constant and < c < 1 as discussed 
previously. The characteristic mode energy w p in Eq. 
( 25 ) is initially chosen as 



Vp 2 



=Q 2 



(27) 



where we added a mass term proportional to the Debye 
screening mass to avoid infrared divergencies. Finally 
we fix the non-abelian Gauss-law constraint by use of 
a standard relaxing algorithm, which we apply to the 
electric field variables. 



[1] U. W. Heinz, AIP Conf. Proc. 739, 163180 (2005). 
[2] J. Berges, J. -P. Blaizot and F. Gelis, |arXiv:1203.2042| 
[hep-ph]. 

[3] K. Dusling, F. Gelis and R. Venugopalan, J. Phys. G G 
38, 124120 (2011); K. Dusling, F. Gelis and R. Venu- 
gopalan, Nucl. Phys. A 872, 161 (2011); F. Gelis and 
R. Venugopalan, Nucl. Phys. A 776, 135 (2006); F. Gelis 
and R. Venugopalan, Nucl. Phys. A 779, 177 (2006). 

[4] P. Romatschke and R. Venugopalan, Phys. Rev. Lett. 
96, 062302 (2006); P. Romatschke and R. Venugopalan, 
Phys. Rev. D 74, 045011 (2006); K. Fukushima and 
F. Gelis, Nucl. Phys. A 874, 108 (2012). J. Berges and 
S. Schlichting, in preparation 

[5] R. A. Janik and R. B. Peschanski, Phys. Rev. D 74, 
046007 (2006); V. Balasubramanian, A. Bernamonti, 
J. de Boer, N. Copland, B. Craps, E. Keski-Vakkuri, 
B. Muller and A. Schafer et al, Phys. Rev. D 84, 026010 
(2011); M. P. Heller, R. A. Janik and P. Witaszczyk, 
Phys. Rev. Lett. 108, 201602 (2012). 

[6] R. Baier, A. H. Mueller, D. Schiff and D. T. Son, Phys. 
Lett. B 502, 51 (2001). 

[7] A. Kurkela and G. D. Moore, JHEP 1112, 044 (2011). 

[8] J. -P. Blaizot, F. Gelis, J. Liao, L. McLerran, R. Venu- 
gopalan, Nucl. Phys. A 873, 68 (2012). 

[9] A. Kurkela, G. D. Moore, JHEP 1111, 120 ( 2011); 
[10] J. Berges, S. Schlichting and D. Sexty, |arXiv:1203.4646 
[hep-ph]. 

[11] J. Berges, S. Scheffler and D. Sexty, Phys. Lett. B 681, 
362 (2009). 

[12] V. Zakharov, V. Lvov, and G. Falkovich, Kolmogorov 
Spectra of Turbulence, Wave Turbulence (Springer- 
Verlag, Berlin Heidelberg New York, 1992); 
U. Frisch, Turbulence: the legacy of A. N. Kol- 
mogorov, Cambridge Univ. Press, Cambridge, 1995; 
S. Nazarenko, Wave Turbulence, Lecture Notes in 
Physics, Springer, Berlin, 2011. 



[13] J. Berges, A. Rothkopf, J. Schmidt, Phys. Rev. Lett. 

101, 041603 (2008). J. Berges and D. Sexty, Phys. Rev. 

D 83, 085004 ( 2011). J. Berges and D. Mesterhazy, 
|arXiv:1204.1489"l [hep-ph]. 
[14] R. Micha and I. I. Tkachev, Phys. Rev. Lett. 90, 121301 

(2003); R. Micha and I. I. Tkachev, Phys. Rev. D 70, 

043538 (2004). 

[15] Y. Kagan, B. V. Svistunov and G. V. Shlyapnikov, Sov. 

Phys. JETP 84, 279 (1992). 
[16] D. V. Semikoz and I. I. Tkachev, Phys. Rev. D 55, 489 

(1997) 

[17] B. Nowak, D. Sexty and T. Gasenzer, Phys. Rev. 
B 84, 020506 (2011); B. Nowak, J. Schole, D. Sexty 
and T. Gasenzer, Phys. Rev. A 85, 043627 (2012); 
M. Schmidt, S. Erne, B. N owak, D. Sexty and 
T. Gasenzer, arXiv: 1203.3651 [cond-mat. quant-gas]; 
J. Schole, B. Nowak and T. Gasenzer, |arXiv:1204.2487| 
[cond-mat. quant-gas]; B. Nowak and T. Gasenzer 
arXiv: 1206.3181 [cond-mat .quant-gas] . 

[18] P. Arnold and G. D. Moore, Phys. Rev D 73. 0205006 
(2006); P. Arnold and G. D. Moore, Phys. Rev D 73, 
025013 (2006). 

[19] A. H. Mueller, A. I. Shoshi and S. M. H. Wong, Nucl. 

Phys. B 760, 145 (2007); 
[20] M. E. Carrington and A. Rebhan, Eur. Phys. J. C 71, 

1787 (2011). 

[21] J. Berges and D. Sexty, Phys. Rev. Lett. 108, 161601 
(2012); 

[22] J. Berges, S. Scheffler, D. Sexty, Phys. Rev. D 77 (2008) 
034504; J. Berges, D. Gelfand, S. Scheffler, D. Sexty, 
Phys. Lett. B 677 (2009) 210. 

[23] A. Cucchieri and T. Mendes, Nucl. Phys. B 471, 263 
(1996). 

[24] J. Berges, K. Boguslavski and S. Schlichting, Phys. Rev. 
D 85, 076005 (2012); K. Dusling, T. E pelbaum, F. Gelis 
and R. Venugopalan, arXiv: 1206.3336 [hep-ph]; Y. Hatta 
and A. Nishiyama, arXiv: 1206.4743 [hep-ph]. 



